# R script for paper Manger - Pickup
# The Coevolution of Trade Agreement Networks and Democracy
# Journal of Conflict Resolutionracy

#setwd("pathname /Replication Files")

# read data
net73 <- as.matrix(read.table("CoevPTA1973.dat"))
net74 <- as.matrix(read.table("CoevPTA1974.dat"))
net75 <- as.matrix(read.table("CoevPTA1975.dat"))
net76 <- as.matrix(read.table("CoevPTA1976.dat"))
net77 <- as.matrix(read.table("CoevPTA1977.dat"))
net78 <- as.matrix(read.table("CoevPTA1978.dat"))
net80 <- as.matrix(read.table("CoevPTA1980.dat"))
net81 <- as.matrix(read.table("CoevPTA1981.dat"))
net83 <- as.matrix(read.table("CoevPTA1983.dat"))

democ  <- as.matrix(read.table("udsInt73-83_1.dat", na.strings = "NA", as.is = FALSE))
democ1<-cbind(democ[1:145,1],democ[1:145,2],democ[1:145,3],democ[1:145,4],democ[1:145,5],democ[1:145,6],democ[1:145,8],democ[1:145,9],democ[1:145,11])
invGDP <- as.matrix(read.table("invGDP73-83.dat", na.strings = "NA", as.is = FALSE))
invGDP1<-cbind(invGDP[1:145,1],invGDP[1:145,2],invGDP[1:145,3],invGDP[1:145,4],invGDP[1:145,5],invGDP[1:145,6],invGDP[1:145,8],invGDP[1:145,9],invGDP[1:145,11])
GDPcap <- as.matrix(read.table("GDPcap73-83.dat", na.strings = "NA", as.is = FALSE))
GDPcap1<-cbind(GDPcap[1:145,1],GDPcap[1:145,2],GDPcap[1:145,3],GDPcap[1:145,4],GDPcap[1:145,5],GDPcap[1:145,6],GDPcap[1:145,8],GDPcap[1:145,9],GDPcap[1:145,11])

dist <- as.matrix(read.table("Coevdist.dat"))

Ally73 <- as.matrix(read.table("CovAlly1973.dat"))
Ally74 <- as.matrix(read.table("CovAlly1974.dat"))
Ally75 <- as.matrix(read.table("CovAlly1975.dat"))
Ally76 <- as.matrix(read.table("CovAlly1976.dat"))
Ally77 <- as.matrix(read.table("CovAlly1977.dat"))
Ally78 <- as.matrix(read.table("CovAlly1978.dat"))
Ally80 <- as.matrix(read.table("CovAlly1980.dat"))
Ally81 <- as.matrix(read.table("CovAlly1981.dat"))
Ally83 <- as.matrix(read.table("CovAlly1983.dat"))
        
Trade73 <- as.matrix(read.table("Cvtrade1973.dat"))
Trade74 <- as.matrix(read.table("Cvtrade1974.dat"))
Trade75 <- as.matrix(read.table("Cvtrade1975.dat"))
Trade76 <- as.matrix(read.table("Cvtrade1976.dat"))
Trade77 <- as.matrix(read.table("Cvtrade1977.dat"))
Trade78 <- as.matrix(read.table("Cvtrade1978.dat"))
Trade80 <- as.matrix(read.table("Cvtrade1980.dat"))
Trade81 <- as.matrix(read.table("Cvtrade1981.dat"))
Trade83 <- as.matrix(read.table("Cvtrade1983.dat"))

# Prepare data for RSiena.

library(RSiena)

# create dependent variable
PTAs <- sienaNet(array(c(net73, net74, net75, net76, net77,
               net78, net80, net81, net83),
                        dim = c(145, 145, 9) ))

#create dependent behavior variable 
Democracy <- sienaNet(democ1, type="behavior")

# create monadic time varying covariates
GDP_inv   <- varCovar(invGDP1)

# create constant dyadic covariate
Distance  <- coDyadCovar(dist)

# create dyadic time varying covariates
Alliance  <- varDyadCovar(array(c(Ally73, Ally74, Ally75, Ally76, Ally77,
               Ally78, Ally80, Ally81),
                          dim = c(145, 145, 8) ))

Trade   <- varDyadCovar(array( c(Trade73, Trade74, Trade75, Trade76, Trade77,
               Trade78, Trade80, Trade81),
                        dim = c(145, 145, 8) ))


# define data                          
PTAbehdata <- sienaDataCreate(PTAs, Democracy, GDP_inv, Distance,
                           Alliance, Trade)
# create effects structure
PTAeff <- getEffects(PTAbehdata )

# give global report about data
print01Report(PTAbehdata, PTAeff, modelname = 'PTA_coev7383_uds_int' )
   

# Define the kind of non-directed model.
model33 <- sienaModelCreate(projname = "PTA_coev7383_uds_int", modelType = 3, n3=9000)

# Specify the model
PTAeff <- setEffect(PTAeff, transTriads, include=FALSE)
PTAeff <- includeEffects(PTAeff, nbrDist2)
PTAeff <- includeEffects(PTAeff, X, interaction1="Distance")
PTAeff <- includeEffects(PTAeff, X, interaction1="Alliance")
PTAeff <- includeEffects(PTAeff, altX, interaction1="Democracy")
PTAeff <- includeEffects(PTAeff, egoXaltX, interaction1="Democracy")
PTAeff <- includeEffects(PTAeff, altX, interaction1="GDP_inv")
PTAeff <- includeEffects(PTAeff, egoXaltX, interaction1="GDP_inv")
PTAeff <- includeEffects(PTAeff, X, interaction1="Trade")
PTAeff <- includeInteraction(PTAeff, X, altX, interaction1=c("Trade", "GDP_inv"))
PTAeff <- includeEffects(PTAeff, name = "Democracy", avAlt, outdeg, interaction1 = "PTAs")
PTAeff <- setEffect(PTAeff, name = "Democracy", quad, include=TRUE)

# Check the current model specification
PTAeff

# estimate
ans1 <- siena07(model33, data = PTAbehdata, effects = PTAeff)
ans1 
